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The infinite Projected Entangled Pair States (iPEPS) algorithm [J. Jordan et al, PRL 101, 250602 
(2008)] has become a usefnl tool in the calculation of ground state properties of 2d quantum lattice 
systems in the thermodynamic limit. Despite its many successful implementations, the method has 
some limitations in its present formulation which hinder its application to some highly-entangled 
systems. The purpose of this paper is to unravel some of these issues, in turn enhancing the stability 
and efficiency of iPEPS methods. For this, we first introduce the fast full update scheme, where 
effective environment and iPEPS tensors are both simultaneously updated (or evolved) throughout 
time. As we shall show, this implies two crucial advantages: (i) dramatic computational savings, 
and (ii) improved overall stability. Besides, we extend the application of the local gauge fixing, 
successfully implemented for finite-size PEPS [M. Lubasch, J. Ignacio Cirac, M.-C. Banuls, PRB 
90, 064425 (2014)], to the iPEPS algorithm. We see that the gauge fixing not only further improves 
the stability of the method, but also accelerates the convergence of the alternating least squares 
sweeping in the (either “full” or “fast full”) tensor update scheme. The improvement in terms of 
computational cost and stability of the resulting “improved” iPEPS algorithm is benchmarked by 
studying the ground state properties of the quantum Heisenberg and transverse-field Ising models 
on an infinite square lattice. 
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I. INTRODUCTION 

In the last years, tensor networks have emerged as a 
powerful tool to understand quantum many-body sys¬ 
tems.® From the point of view of numerical simulations 
there have been a number of novel algorithms developed, 
whose inner workings are deeply rooted in the theory of 
quantum entanglement. In all these algorithms quantum 
many-body states are conveniently represented by ten¬ 
sor networks that efficiently capture the natural struc¬ 
ture of quantum correlations in the swtem (such as, e.g., 
the so-called entanglement area /aw.l^HSl). Many impor¬ 
tant classes of states can be accurately approximated by 
a tensor network with a number of parameters that de¬ 
pends only polynomially on the size of the system. Such 
properties, among others, have enabled tensor network 
methods to break the curse of dimensionality, namely, 
the fact that the Hilbert space dimension of a quantum 
many-body system is exponentially large in the number 
of particles. Thanks to this, it is now possible to effi¬ 
ciently simulate quantum many-body systems by target¬ 
ing the relevant tiny corner of quantum states (e.g., those 
satisfying an area-law) inside of the exponentially-large 
Hilbert space.®® 

The so-called matrix product state (MPS]p®^ is a 
typical tensor network ansatz representing the state of 
Id gapped quantum lattice systems. It is also well 
known that MPS is the class of variational wave func¬ 
tions at the root of the de nsity matrix renormaliza¬ 
tion group (DMRG) method, ^21121 widely used in the 
study of Id systems.lSH^ Subsequently, MPS methods 


to study time-evolution of Id systems have also been 
put for ward , such as time evolvin g blo ck decimation 
(TEBD),®^ time-dependent DMRGp®^ and more re¬ 
cently algorithms based on the time-dependent varia¬ 
tional principle.!^ One way of generalizing MPS methods 
to higher dimensional systemJ^ is using projected en¬ 
tangled pair states (PEPS),sometimes also called 
tensor product states (TPS).l22H3l] There has been a lot 
of progress both in conceptual and algorithmic develop¬ 
ments for PEPS, see, e.g., Refs [521 - H51 From the numer¬ 
ical perspective, PEPS have been used to study ground 
state properties as well as dynamics of 2d lattice systems, 
both of finite and infinite size. Moreover, motivated by 
the success of Id m ethods in the thermodynamic limit 
such as iTEBD, 1^2150] so-called infinite-PEPS (iPEPS) 
algor it hirPSmi was put forward to study infinite-size 2d 
quantum lattice systems. 

So far, the iPEPS algorithm has been quite successful 
in studying ground state properties of a growing number 
of 2d quantum lattice systems (see, e.g., Ref. [T]and ref¬ 
erences therein). In general terms, results obtained by 
using iPEPS can be competitive when c ompa red to the 
ones derived from quantum Monte Carlo.l^iEl] And what 
is more important, the iPEPS algorithm is not hampered 
by the sign problem (unlike quantum Monte Garlo) when 
studying fermionic and frustrated spin systems. Recent 
applications of iPEPS to such systems include calcula- 
tions for the t — J model of fermions on the square ,1^3111111 
and honeycomb lattices,!^ as well as the J 1 — J 2 frustrated 
Heisenberg model on the square lattice,^ the Shastry- 
Sutherland model,^ and the Kagome Heisenberg anti- 
ferromagnet .1^ 
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One of the main drawbacks of the iPEPS algorithm 
is its high computational cost as a function of the 
bond dimension D which controls the accuracy of the 
method. This is particularly true when trying to ob¬ 
tain good accuracies in physical regimes where entan¬ 
glement is large (e.g., close to a quantum critical point, 
or in the presence of many nearly-degenerate quan¬ 
tum states), which requires a large D. The compu¬ 
tational bottleneck of the method is the calculation of 
the so-called effective environment, i.e., the effective 
description of the tensor network surrounding a given 
site. Technically, effective environments can be com¬ 
puted using different approach es, such as TRG/SRG 
and HOTRG/HOSRG^^n^ iTEBD,!^ corner trans¬ 
fer matrices (GTM),I^3II1H63] qj. recently the ten¬ 

sor network renormalization method.!^ Independently of 
the chosen approach for their calculation, accurate ef¬ 
fective environments are required for the so-called full 
update (FU), which is the accurate scheme proposed to 
find the iPEPS tensors throughout a (imaginary) time 
evolution. Because of the large computational cost only 
relatively small bond dimensions can be afforded when 
using the FU scheme. 

Alternative update schemes have been developed try¬ 
ing to overcome this problem, but only with partial suc¬ 
cess. A ver y pop ular approach is the so-called simple 
update (SU),I25E21 which relies on a mean-held approxi¬ 
mation of the effective environment, thus being very ef- 
hcient. Such a scheme allows to reach large bond di¬ 
mensions in the iPEPS but, not surprisingly, it does not 
produce accurate results when systems are very strongly- 
correlated. Intermediate approaches interpolating be¬ 
tween the FU and the S U sch emes have also been put 

forward as an alternativel '^'^FS I 

So, here is the dilemma: accurate update schemes like 
the FU are too costly, whereas efficient update schemes 
such as the SU are not accurate enough. The question 
then is: can one somehow “accelerate” the FU, making 
it more efficient while keeping its accuracy? 

In this paper we give a positive answer to this question. 
We do this by constructing an update scheme, which we 
call “fast full update” (FFU), that signihcantly reduces 
the computational cost of iPEPS algorithms while still 
being accurate. More specifically, in this new strategy at 
every time step the tensors of the effective environment 
are updated by a single iteration step (in a sense to be 
made specific later) and simultaneously with those of the 
iPEPS. Importantly, we find that applying this strategy 
to iPEPS algorithms not only reduces the computational 
cost by a large factor (as expected), but also contributes 
to stabilize the algorithm. The reason for this is that the 
successively updated environment helps to maintain the 
compatibility between the related tensors throughout the 
time evolution, as we shall explain later. 

In addition, we show that incorporating the local gauge 
fixing scheme proposed in Ref. |46| (and successfully ap¬ 
plied to finite PEPS) can make the iPEPS algorithm even 
faster and more stable. As described in Ref. |3S1 the idea 


for the gauge fixing of PEPS tensors is inspired by the 
case of MPS, for which tensors can always be represented 
in a canonical form during their update by means of lo¬ 
cal gauge fixing. In the canonical form many of the ten¬ 
sor manipulations of an MPS get simplified (or directly 
canceled out) implying a much better conditioning and 
stability of related algorithms. However, unlike in the 
MPS case, there is no exact canonical form for PEPS in 
the same sense. A recent attempt along thi s directi on is 
the so-called quasi-canonical form for iPEPS.^^^EUlZlxj^ig 
has been shown to lead to some computational advan¬ 
tages, but unfortunately does not fully capture the effect 
of quantum correlations spreading throughout the 2d lat¬ 
tice. A different approach was considered in Refs. 1461471 
where it was shown that by considering the effect of the 
entire 2d lattice, a local gauge choice of the tensors can 
also produce a well-conditioned environment, which in 
turn improves the stability of the subsequent calcula¬ 
tions. Here we apply the same local gauge fixing as in 
Ref. Uni to the iPEPS algorithm and show that it not 
only improves the stability, but also accelerates the con¬ 
vergence of the alternating least squares sweeping in the 
tensor update scheme. 

To show the validity of these approaches, we provide 
benchmarking calculations for the “improved” iPEPS al¬ 
gorithm with the two “improvements” mentioned above 
(FFU -I- gauge fixing). In particular, we analyse the 
computational cost and the stability of the algorithm, 
for ground-state calculations of the Heisenberg and 
transverse-field Ising models on an infinite square lat¬ 
tice. We shall see quantitatively that the “improve¬ 
ments” both accelerate and stabilize the overall numeri¬ 
cal calculations. 

The paper is structured as follows. Some background 
material on PEPS, iPEPS, and the iPEPS algorithm is 
presented in Sec. H. In Sec. HI we introduce the im¬ 
provements (FFU and gauge-fixing). Benchmarking cal¬ 
culations are presented in Sec. IV. Finally, conclusions 
are presented in Sec. V. 

II. BACKGROUND 

A. PEPS and iPEPS 

1. Generalities 

For completeness, we briefly review the notation and 
fundamental properties of projected entangled pair states 
(PEPS).l2SH2ll To this end, let us consider a 2d quantum 
lattice system consisting of N sites, each of which is de¬ 
scribed by a local Hilbert space C"^. The full Hilbert space 
of the system is thus H = (C*^)®-^. As the dimension of 
the full Hilbert space grows exponentially with the size 
of the system, the problem quickly becomes intractable 
already for moderately-low values of N. In order to avoid 
this curse of dimensionality, an option is to use a PEPS to 
represent a pure state of the system. Generally speaking. 













a PEPS is a state defined by a 2d lattice of interconnected 
tensors, i.e., 
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where is the local basis of the site located at jy = 
{xi,yi). Depending on the geometry of the lattice pat- 
tern, the tensor As~^ at each lattice site contains nf*. 
bond indexes taking up to D values (nf. is typically the 
number of nearest neighbours of the lattice site ri) and 
a physical index taking up to d values. The operation 
F contracts all the tensors along the bond indices. 
Conventionally, D is usually referred to as the bond di¬ 
mension, which plays the role of a parameter quantifying 
both the size of the tensors in the PEPS, and also the 
amount of entanglement in the wave function.l^ Thus, 
the larger D, the better the PEPS can represent the state 
of the system (since there are more variational parame¬ 
ters). As an example, in Fig. [^a) we illustrate a graphi¬ 
cal representation of the PEPS for a 5 x 5 square lattice. 
In this case the amount of complex parameters describing 
the PEPS quantum state is 0{NdD'^), thus polynomial 
in N, D and d, in contrast with the exponential depen¬ 
dence on N for an arbitrary state in the Hilbert space. 
For the sake of simplicity, from now on we shall always 
assume that we have a 2d square lattice.!^ 


2. Properties 

PEPS have some important properties that make it an 
appropriate representation for 2d quantum lattice sys¬ 
tems. First, and similarly to MPS, PEPS satisfy the 
entanglement area law.^HS More specifically, the scaling 
of entanglement entropy of an L x L block of a PEPS 
is 0{L\ogD). Accordingly, PEPS can describe very well 
the entanglement structure of many interesting 2d quan¬ 
tum systems, including low-energy eigenstates of many 
2d Hamiltonians with local interactions. Second, PEPS 
can in principle represent systems with both finite and 
infinite correlation length.^SlThis is to be contrasted with 
the case of MPS, where only a finite correlation length 
is possible. Third, and also unlike for MPS, given the 
loops present on 2d lattices there is no obvious canonical 
form for a PEPS (although some prop osals ha ve been re¬ 
cently put forward along this directiorP^IllIlZl) _ Last, but 
not least, PEPS can be used to represent systems in the 
thermodynamic limit by using a small number of tensors, 
under the assumption of shift invariance. More precisely, 
a unit cell of tensors is repeated all over the 2d lattice to 
construct an arbitrarily-large shift-invariant PEPS. For 
an infinite system, this is the so-called infinite-PEPS, or 
iPEPS (see Fig. Bb) for a graphical illustration of an 
iPEPS with a two-site unit cell). 



Figure 1: (Color online) (a) Graphical representation of a 
PEPS on a 5 X 5 square lattice. Each ball is a tensor, and 
lines correspond to tensor indices. Lines from one tensor to 
another correspond to common summed indices (or contracted 
indices). The free index, or open leg, in each tensor is called 
physical index, and corresponds to the local degrees of freedom 
of the local Hilbert space at every site, (b) An infinite PEPS 
with a two-tensor unit cell. The two tensors A and B are 
repeated on the infinite 2d lattice. 


3. Numerical application 

PEPS can be used to study both ground state proper¬ 
ties as well as dynamics of 2d quantum lattice systems. 
For ground states, one can either (i) variationally opti¬ 
mize the PEPS tensors so as to minimise the expectation 
value of the 2d corresponding Hamiltonian (as done in 
DMRG in Id), or (ii) evolve the system in imaginary time 
until a fixed point (ground state) is reached (as done in 
TEBD in Id). The second method can also be applied to 
study real time evolution of the system. This approach 
is also easy to extend to th e thermodynamic limit, called 
the iPEPS algorithmwhich we go over again briefly 
in the next section. 


B. The iPEPS algorithm 

1. Generalities 


For a given Hamiltonian H, the ground state of the 
system can be obtained by evolving an initial state 40) 
in imaginary time j3 as described by 


4gs) 


lim 
/3—>-oo 


e-^^4o) 

e-^h\^o)\\' 


( 2 ) 


Moreover, the real-time evolution is described via the 
solution of the Schrodinger equation, which for a time- 
independent Hamiltonian FI reads 


4(t)) = e-*"‘4o). 


( 3 ) 
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From now on, let us consider the evolution of an iPEPS 
in imaginary time (the extension to real-time evolution 
is technically straightforward). 

We assume that the Hamiltonian contains only trans- 
lationally invariant nearest-neighbour interactions, i.e., 

iJ = ^ (4) 

{r,r') 

where the sum is performed over all the nearest neigh¬ 
bours {r,r'). For a Hamiltonian that is invariant under 
translations, we can use an iPEPS also with translation 
symmetry to represent the wave function |4') of the sys¬ 
tem. This could be achieved by repeating the same tensor 
on each lattice site all over the lattice. However, our up¬ 
date scheme (described below) requires neighboring ten¬ 
sors to be different, which is why we typically use a two- 
site translationally invariant iPEPS shown in Fig. [^b), 
depending only on two tensors A and B. If translational 
symmetries are spontaneously broken in the thermody¬ 
namic limit, then a larger unit cell of tensors (compatible 
with the structure of the ground state) can be used in¬ 
stead.!^ 

In order to compute the ground state of the system by 
an evolution in imaginary time, one first decomposes the 
time evolution operator into a product of so-called two- 
body gates. To this end, the Hamiltonian H is rewritten 
as 

H = Hi + Hr+ H^+ Hd, (5) 

where each term Hi = X](r f')6i ^ is 

the sum of mutually-commuting Hamiltonian terms for 
links labelled as {left, right, up, down). Notice, though, 
that the commutator [Hi,Hj] is in general different from 
zero whenever i j. Applying the first-order Suzuki- 
Trotter decomposition,!^ we can write the time evolution 
operator as 

« ( 6 ) 

where 5 is an infinitesimal time step and m = /3/(5 G N 
is the number of “time steps” that needs to evolve the 
system, so as to reach total evolution time j3. Since each 
term Hi is a sum of mutually-commuting terms, we can 
further write each term e~^'^ in Eq. (j^ as a product of 
two-body gates, i.e., 

= n (7) 

{r,r')£i 

with 

For a given time step of the evolution, let us consider 
the action of the term e~^'^ on an infinite PEPS |'I'. 4 b) 
with two tensors A and B, and bond dimension D. To 
this aim we focus on one link on the lattice and, only 
initially, disregard the effect of the gates gt’"’’" 1 on the 


rest of the links (which is approximately correct, since 
(5 <C 1 and thus g!”’” 1^1). Let us consider the case of 
an r—link, for concreteness. After applying the gate, we 
obtain a new iPEPS \^a'B') which is characterized by 
tensors A and B everywhere except for the two tensors 
connected by the link where the gate acted. More pre¬ 
cisely, 14'^/^/) and |4 '^b) differ from each other only by 
two tensors. Because of the effect of the gate, the bond 
dimension of the affected index changes to D' < d^D, 
and thus increases, corresponding to a change of the en¬ 
tanglement in the tensor network. Because of this, the 
iPEPS bond dimension quickly increases exponentially 
fast after few gate applications, making the simulation 
intractable. To overcome this problem, the infinite PEPS 
14'^/^/) is approximated by a new PEPS I'l'^jj), by re¬ 
placing tensors A! and B' by A and B, where these two 
last tensors have again bond dimension D for the affected 
index. This is done in a way such that the state |'I'^^) 
is close to the exact state 14'^/^/), and thus introduces a 
small error only. Such a procedure is called tensor update 
of the iPEPS. 

A possibility to implement the tensor update is to look 
for new tensors A and B that minimise the squared dis¬ 
tance between the exact and the approximating state, 

i.e. 

min|||4'A'B') - |4 'ab)IP = niind(i,i?), (8) 

A,B A,B 

with 

d{A,B) = (4-^,b,|«'a'B') + (4'^bI%b)- 

-{^abI^A'b') - {'^abI^A'b')- (9) 

To solve the problem in Eq. (|^, two main tasks are 
needed. These are (i) the effective environment calcu¬ 
lation, and (ii) the tensor update. Since these pieces of 
the method will be fundamental for the improvements 
to be later explained, we review them in detail in the 
following. 


2. Effective environment calculation 

To properly evaluate d{A, B) one needs to take into 
account the effect of the whole tensor network surround¬ 
ing the affected link, i.e., the environment. Such a ten¬ 
sor network of infinitely-many tensors is conveniently ap¬ 
proximated by an effective environment, consisting of 
a small number of tensors only. The effective environ¬ 
ment can be computed using v arious app roaches, such as 
TRG/SRG, H0TRG/H0SRG,I22MI1IM1 tensor network 
renormalization (TNR),!^ iTEBD,!^ and corner transfer 
matrix (GTM) methods. This last approach will be our 
choice in this paper. We shall not explain GTM methods 
in full detail here, and we refer the reader to the extensive 
existing literature on the topic such as Refs. P1M I55] 
for technicalities. However, since these methods will also 
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Figure 3: (Color online) Tensor network diagrams showing 
how to compute T, R and S in Eq. (12 I. 


3. iPEPS Tensor Full Update 

Eq. ^ implies that the squared distance d{A,B) is 
a quadratic function of tensors A and B. Thus, once 
the effective environment is obtained, one can apply an 
alternating least squares method to optimize tensors A 
and B. This is the so-called full update (FU), and it 
follows the steps below: 

(i) Fix tensor B to some initial tensoil^or to the ten¬ 
sor obtained from previous iteration. In order to 
find A, rewrite Eq. ([^ as a quadratic scalar ex¬ 
pression for the tensor, 

= A^RA-A'fS - S^A + T. (12) 


Figure 2: (Color online) (a) Left: 2d lattice of tensors formed 
from a and 6; right: contractions to obtain tensors a and b. 
(b) Environment of a given link on the lattice (here an r- 
link). (c) Effective environment of a given link on the lattice 
(here an r-link). (d) 6-tensor representation of the effective 
environment around the link. 


be important at a later stage of this paper, we review 
briefly some notations and conventions. 

We first construct an infinite square lattice C{a,b) by 
contracting the physical indexes of and (iPyiBl, see 

Fig. I^a) where 


^Irud 

^Irud 


^Irudi^l'r'u'd’T 

s 

( 10 ) 

Blrudi^l'r'u'd')* J 

( 11 ) 


s 


where I, r, u, d are combined indices, e.g. I = (/, I'). 

The exact environment £{ri,r 2 ) of two sites at 
r*! and r 2 is shown in Fig. I^b). In CTM meth¬ 
ods, this is approximated by an effective tensor net¬ 
work G{ri,r 2 ), the effective environment, which com¬ 
prises a set of four x x y corner transfer matrices 
{Cl, (72, Cs, ( 74 }, eight half transfer row/column tensors 
{Tai,Ta 2 ,Ta 3 ,Ta 4 ,Tbi,Tb 2 ,Tb 3 ,Tb 4 } and two tensors 
a and b, see Fig. Sc). A further simple contraction 
of the tensor network produces an effective environ¬ 
ment for the considered link in terms of only six tensors 
{Ei,E2, E3, E4, F 5 , Eq}, see Fig. [^d). 


In this equation, A is understood as a reshaped vec¬ 
tor with D'^d components and matrices R, S, T can 
be obtained from the appropriate tensor contrac¬ 
tions including the effective environment around 
the r-link, see Fig. 


(ii) We find the minimum of d{A, A^) in Eq. ( 12 ), with 
respect to A\ which is given by A = R~^S. 


(hi) Next, we fix the tensor A and find B using the same 
procedure as steps (i) and (ii) above. 


The above steps are iterated until the cost function 
d{A, B) converges to a sufficiently-small value. A possi¬ 
bility to check convergence is, e.g., to check the value of 
this cost function between two successive iterations and 
compare it to some small tolerance. Last but not least: 
once the optimal tensors are found, these are replaced 
over the entire 2d lattice approximating the effect of all 
the gates ^ acting over the infinitely-many links of 
the same type. Such a procedure defines the updated 
infinite PEPS I'L^s) terms of the two new tensors. 
Finally, the same procedure is repeated for the 1-, u- and 
d-links to complete one time step. This is iterated un¬ 
til the desired real running time is achieved, or until the 
iPEPS converges to a fixed-point approximation of the 
ground state for imaginary-time evolutions. 

As such, the computational cost of the procedure above 
is quite high due to the calculation of the inverse of the 
matrix R in step (ii). Specifically, this cost is 0{D^) if 
recurrent methods are applied for computing the inverse 
(e.g., biconjugate gradient). Otherwise, the cost would 
be 0{D^^) for an exact inverse calculation. 
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Figure 4: (Color online) (a) Subtensors for iPEPS tensor A 
and B. Decompose A = Xaj{ and B = by means of 
decompositions such as QR or SVD. (b) The action of the 
two-body gate g on the iPEPS tensors A and B connected by 
an r-link is equivalent to its action on the subtensors an and 
bn only, leaving X and Y unaffected. 


4 . Reduced Tensor Full Update 


As an efficient and convenient alternative to this 
method, one can apply the “revised” FU scheme dis¬ 
cussed in Ref. [37] where, instead of updating two ten¬ 
sors A and B, one updates some lower-rank subtensors 
related to them (sometimes they are also called reduced 
tensor^^. These subtensors are denoted as qr and bn, 
respectively. The main idea of this optimization scheme 
is based on the observation that applying the two-body 
gate g on two sites A and B connected by a specific link 
changes the properties of that link only^ while the oth¬ 
ers remain unchanged. For instance, for an r-link this 
results into modifying the bond dimension for this link 
from D to D' > D, while the size of the other indices is 
unaffected. 

Thus, by means of the QR or singular value decompo¬ 
sition (SVD), we decompose tensors A and B such that 
A = Xor and B = bpY, where qr and bp are connected 
by the bond index on the r-link, and also contain the 
physical indices (see Fig.j^a)). When applying the two- 
body gate g, we now contract it with or and bp in order 
to update them directly. Once qr and bp are updated 
to, say, clr and bn, we can easily get updated iPEPS ten¬ 
sors A and B by using A = Xor and B — bnY. The 
update scheme for or and bn can be performed similarly 
as we explained for A and B, following steps (i-iii) above. 
More precisely, in these steps we replace A and B by cir 
and &L, which are alternatively optimized according to a 
similar figure of merit defined as 

dCanM) = + 

Eq. ( |12| ) for ur is now rewritten as, 

d{dR, oj^) = d^Rd_R - - S^or + T, (14) 


and the cost function for variable bp is defined in a sim¬ 
ilar way. Note that the tensors {R,S,T} in Eq. 14 are 
different from the ones in Eq. see Eig. 

The computational cost of this update scheme is 
0{d^D^), where the inverse of the matrix R can now be 



Figure 5: (Color online) Tensor network diagrams showing 
how to compute T, R and S in Eq. (141. 


computed exactly. Due to this huge advantage in com¬ 
putational cost, as compared to the direct update of the 
iPEPS tensors A and B, the present scheme is able to 
deal with iPEPS of larger bond dimension D, with the 
corresponding advantages. 

The rest of the update follows as explained in the pre¬ 
vious subsection. This is, the reduced tensors are op¬ 
timized until the cost function is sufficiently small, and 
then, these are replaced over the entire 2d lattice. As 
explained before, this last step approximates the effect 
of all the gates acting on all the links of the lattice of the 
same type. As before, the procedure defines an updated 
infinite PEPS in terms of two new tensors. After 

this, one repeats the same procedure for the 1-, u- and 
d-links to complete one time step of the evolution. This 
will be iterated until the iPEPS converges to a fixed point 
approximating the ground state (for imaginary time), or 
until the desired real-time evolution is completed. 


III. IMPROVING THE IPEPS ALGORITHM 

One of the main limitations of the FU iPEPS algo¬ 
rithm in its present formulation is that the effective en¬ 
vironment has to be recomputed from scratch at ev¬ 
ery time step. Obtaining a converged effective environ¬ 
ment requires Actm iterations of the CTM algorithm 
with Nctm depending on the amount of entanglement 
in the system (for not-too-entangled systems, typically 
VcTM ~ 10 — 20, but it can be considerably larger in 
strongly entangled systems). Since these CTM iterations 
are the computational bottleneck of the method, it is de¬ 
sirable to keep VcTM as low as possible. However, if one 
uses a VcTM which is too small, then the errors intro¬ 
duced in the effective environment can lead to instability 
problems of the method. 

Here we propose two improvements to the algorithm, 
which enhance the stability of the method and make it 
more efficient. Eirst, we explain how to implement what 
we call a fast full update (FFU), where the accuracy of the 
EU is preserved while substantially reducing its computa¬ 
tional cost and improving its stability. Second, we discuss 
the application of a local gauge fixing, as in Ref. 1461 which 
naturally improves stability and also accelerates the con¬ 
vergence of the overall method. This is explained in what 
follows, and benchmarking numerical calculations shall 
be provided in the forthcoming sections. 
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A. Fast Full Update (FFU) 

In this update scheme we reduce the computational 
cost by using the following idea: instead of recomput¬ 
ing the environment tensors from scratch we can update 
the environment tensors simultaneously with those of the 
iPEPS at each step, just by a single CTM iteration step. 
This reduces the computational cost by a large factor 
(see results below). The crucial technical point of this 
FFU scheme is to make sure that the environment ten¬ 
sors remain compatible with the updated iPEPS tensors 
throughout the (imaginary) time evolution. 

The details of the FFU scheme for iPEPS tensors A 
and B plus the effective environment tensors are shown in 
Fig.§ In the following we explain how to proceed at each 
time step, when the two-body gates g are successively 
applied on the r-, /-, u- and d-links: 

(a) r-link update: Suppose that the effective environ¬ 

ment around four links of the iPEPS is charac¬ 
terized by tensors see Fig. |^a- 

i). Applying the gate g to the r-link will mod¬ 
ify the properties of this link, so we need to up¬ 
date the tensors A and B using the update scheme 
from the previous section. Specifically, the effec¬ 
tive environment of an r-link is obtained by first 
absorbing a row of tensors {ra 4 ,5, a, T 62 } to the 
bottom edge and then contracting the tensors ap¬ 
propriately, such that the effective environment is 
represented by six tensors as in Fig.|^a-ii). 

We then apply the tensor update scheme explained 
previously in order to hnd some updated tensors 
Ai and Bi. 

In order to prepare for the next update of the 
link, we now update the environment tensors. 
To this end, we insert two columns of tensors 
{Toi, 5i, oi, T 63 } and {Tbi,ai, 61 , Toa} in the mid¬ 
dle of the tensor network, see Fig. |^a-iii,iv). Note 
that the two tensors a and b connected by link f are 
now replaced by the updated ones oi and bi. But 
importantly, there are still two sites where we use 
the old tensors a and 6 , since the corresponding r- 
link is connected to the old environment, and thus 
was not “formally” updated. 

Next, the columns of tensors {Toi, 61 , a, and 
{r5i, oi, 6 , Toa} are absorbed to the right and left 
edges respectively, see Fig. [^a-iv). After this ab¬ 
sorption, one finds appropriate isometries to renor¬ 
malize the environment tensors as usual, in such a 
way that their bond dimensions do not grow.l^The 
new environment tensors of the /-link are denoted 
as in Fig. |^b-i). 

(b) l-link update: For the action of a gate on the /-link, 
we proceed in the same way as with the r-link, 
which is shown graphically in Fig. ib). The iPEPS 
is now updated and represented by two new tensors 


(a), r-link update 



(h). l-link update 



(c). u-link update 



(d). d-link update 




Figure 6: (Color online) Tensor network diagrams for the 
Fast Full Update. Details are explained in the main text. 


A 2 and B 2 . The new environment tensors for up¬ 
dating the u-link are denoted as in the Fig.[^c-i). 

(c) u-link update: From the tensor network in Fig.j^c- 
i), we obtain the environment tensors for the u- 
link as shown in Fig. [^c-ii). An update of the 
tensors will now produce two new tensors A^ and 
B 3 . We then compute the effective environment for 











































the update of the d-link as shown in Fig. |^c-iv) 

(d) d-link update: Now we follow the same procedure 
as for the u-link in order to update the iPEPS ten¬ 
sors as well as the effective environment. This is 
shown in Fig. id). Finally, we obtain the new 
iPEPS tensors A 4 and B 4 and the new effective en¬ 
vironment of the four links is represented by tensors 
{Cr,Ta^,Tbnti,see Fig.^d-v). 

The above update scheme is illustrated for one time 
step only. One then needs to iterate this scheme in order 
to evolve the system up to the desired time. As in the 
update approaches explained in the previous section, for 
the FFU is also a good idea to choose properly the initial 
state in ground-state calculations (e.g. converged iPEPS 
obtained from the SU scheme, or from the FU scheme 
with a smaller bond dimension). Evidently, for the FFU 
this choice also helps in the stability and fast convergence 
of the algorithm. 

The FFU that we just presented has two key advan¬ 
tages: first, one keeps an environment at every step that 
is perfectly compatible with the tensors in the iPEPS in 
all bond indices. This is the reason why, e.g., in Fig.|^a- 
iv) we still have two tensors a and b at two sites. Such 
a property naturally improves the stability. Second, the 
environment is not reconverged for every link at every 
stepV^ As a result, this reduces the required computa¬ 
tional time considerably. 


B. Gauge fixing 


In contrast to MPS, a PEPS does not have a canon¬ 
ical form. Therefore, it is difficult to fix the gauge of 
the PEPS tensors in some appropriate way during the 
time evolution. Despite this, local gauge fixing schemes 
have pr oven usefu l in improving the stability of the al¬ 
gorithm In this paper, we use the local gauge 

fixing proposed in Ref. HU (applied there to finite PEPS), 
in the context of the iPEPS algorithm. We shall see that 
this not only helps to improve the stability of the method, 
but also to accelerate its convergence. 

To this end, we consider the tensor update applied to 
the reduced tensors. We rewrite Eq. (13) as follows: 


diflR, II) = a'lb'lNLRaRbf + d}^\NLRdRbR - 

-^Wh^LRa'Rbf - a'j^LAfLRa'Rbf.ilS) 


where Mlr is the “norm” tensor obtained by comput¬ 
ing the overlap of l^ajj&j,) '''^hile leaving out the 

tensors or, b^ and their complex conjugates. Specifically, 
this can be done by contracting the tensor network shown 
in Fig. (a). Also, the cost function is represented di- 
agrammatically in Fig. (b). In order to update the 
subtensors one needs, thus, to compute the norm ten¬ 
sor AfLR. 

Tensor N'lr has the following properties: first, it is im¬ 
possible to choose a gauge in such a way that this tensor 



Figure 7: (Color online) (a) Contraction producing the norm 
tensor. The leading computational cost of this contraction is 
0{d*D*X^) + 0{d^D^x^)+Oid^D^x^f^- (b) Diagrammatic 
representation of the cost function defined in Eq. (151. 



Figure 8: (Color online) (a) Positive approximant for the 
norm tensor Nlr- (b) Apply the QR decomposition for the 
left and the right bond indices of tensor Z, so that Z = 
QlR = LQr. (c) Insert the identities 7 = L~^L = R~^R 
into the left and right bond indexes of the tensor Z to obtain 
Z = L~^ZR~^ and the new subtensors dn = Lur, 6 _l = 6 l 7 ?. 
(d) A new norm tensor is obtained as JVlr = ZZK In or¬ 
der to keep the compatibility once the subtensors have been 
updated, we need to multiply the tensors X and Y by ma¬ 
trices L~^,R~^, so that X = XL~^ and Y = R~^Y, respec¬ 
tively. This is done before recovering the tensors A = Xclr 
and B = ftnU. 


is the identity matrix for all links at the same time. The 
reason is that an iPEPS does not have a canonical form 
in the same sense as MPS. Second, this tensor needs to 
be computed using approximations, which implies that, 
generally, it is neither strictly Hermitian nor positively 
defined. Although one could always apply some approxi¬ 
mation methods that preserve positivity explicitly (such 
as the single layer method,!^, the tensors obtained from 
such approaches do not produce as accurate results as 
other methods that do not enforce positivity (such as 
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Figure 9: (Color online) Mean value of the relative change 
Ed = Mti+i — rfu|/|rfinit| of the cost function d and the local 
fidelity Sab = - /ai,l/l/™''l in the iPEPS algorithm with 

gauge (filled symbols) or without gauge (open symbols), for 
the Ising model with transverse held (a,b) and the Heisenberg 
model (c,d) respectively. For the Ising model the magnetic 
held is h = 3.01 and {D, x) = (3,30), whereas for the Heisen¬ 
berg model {D,x) = (5,50). The time step is 5 = 0.005. The 
mean value is taken over all steps in imaginary time. 


CTM methods).!^ It is well-known that the (normally 
small) negative part of the norm tensor often causes some 
ill-posed conditions in updating the iPEPS tensors. To 
circumvent this problem, one can make it Hermitian and 
positive-defined using the following stepsP^ first, one ap¬ 
proximates the norm tensor as J^lr = (J^lr +-^Li?)/2, 
which is Hermitian. Next, we approximate Mlr by its 
positive part. To achieve this, one applies the eigen¬ 
value decomposition Nlr = and replaces the 

(small) negative eigenvalues in E by zero. This is the 
so-called positive approximant. The approximate eigen¬ 
values are now denoted Moreover, Nlr = ZZ\ 

where Z = see Fig.|^a). 

The gauge-fixing that we apply here is explained in 
Fig.m After fixing the local gauge in the norm tensor, 
we replace Nlr and compatible subtensors into Eq. (151, 
and then start the variational update of subtensors for 
the iPEPS. As shown in Ref. HHl this choice of gauge im¬ 
proves a lot the conditioning of the norm tensor, and thus 
greatly increases the stability of the update. In partic¬ 
ular, when this gauge-fixing is combined with the FEU, 
the resulting iPEPS algorithm is remarkably fast, stable 
and accurate. 


IV. RESULTS 


We have benchmarked the improved iPEPS algorithm, 
including both the FEU and the gauge-fixing, by study¬ 
ing ground state properties of two models on an infinite 


2d square lattice. The first model is the ferromagnetic 
quantum Ising model with transverse magnetic field, 

H = (16) 

r 

where is the i = {x, z) Pauli matrix for site r, h is the 
transverse magnetic held, and (r, r') represent nearest- 
neighbor sites. The second example is the spin-1/2 anti¬ 
ferromagnetic Heisenberg model, 

(17) 

(r,r') 


where 51^ = (crif^, crlf, tTly')/2. 

In our simulations, we have represented the ground 
state of the system by a two-site translationally invari¬ 
ant iPEPS made up of two tensors, A and B. In or¬ 
der to approximate the ground state of the system, we 
have applied imaginary-time evolution together with a 
second-order Suzuki-Trotter decomposition using 5 down 
to 10“^ — 10“"*. To update the tensors, we used the 
alternating-least-square (ALS) sweep for the subtensors, 
as explained in Sec. |HB4[ combined with the FEU. At 
every update, we have also hxed the gauge of the tensors 
according to the gauge-hxing described above. Let us 
stress here that the leading order of the computational 
cost is the same for both FU and FEU -I- gauge-hxing 
schemes, but the prefactor and the subleading correction 
are different. These turn out to produce a big difference 
in practical running times, as we shall see. 

In order to assess the advantage of the local gauge Hx- 
ing for the norm matrix Nlr. used in the ALS sweep, 
we Hrst compare the mean condition number of this ma¬ 
trix between using and not using the gauge. As a rule 
of thumb, the larger the condition number of Nlr, the 
less accuracy we get in solving the system of linear equa¬ 
tions at every step of the ALS. The result is shown in 
Table. Ufor different models and bond dimensions. Over¬ 
all we see that the condition number of the norm matrix 
in the case of gauge hxing is improved by several orders 
of magnitude when compared to the case without gauge 
fixing. For completeness, we also shown in the table the 
condition numbers of matrices L and R. The gauge fix¬ 
ing, thus, improves the stability of the iPEPS algorithm. 
This result is very similar to what has been obtained for 
finite PEPS in a similar context.l^ 

Besides, we observe that the gauge fixing in the norm 
matrix also accelerates the convergence in the ALS 
sweeping. More concretely, in Fig. we show the conver¬ 
gence of the relative change of the cost function defined 
in Eq. (15), as well as the local fidelity of the subtensors 
between two iterations u and u -I- 1 defined as 


fU + l _ 

ab 


(a 


u+l^^u+l 


albl) 


U+l? R + 1 I 2i+l 7 lA + 1 
Xd Ut \ar> Ut 


){a'ib 


U Uu 

r'^l\^r'^l 


bl) 


(18) 


where, e.g., |a^5^), is to be understood as tensors and 
6 ^ with their indices reshaped as a vector. We observe 
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(Ax) 

No Gauge 



Gauge 


R 




L 




Heisenberg model 

(2,20) 

(9.33 ± 1.50) 

X 

10^ 

1.46 ±0.01 


5.49 ±0.23 



5.29 ±0.39 



(3,30) 

(1.6 ±0.17) X 10® 

(0.15 ±0.00) X 

10® 

(0.41 ± 

0.01) 

X 

10® 

(0.45 ± 

0.01) 

X 

10® 

(4,40) 

(8.36 ±0.30) 

X 

10^ 

(3.02 ±0.02) X 

10" 

(1.78 ± 

0.00) 

X 

10" 

(1.80 ± 

0.02) 

X 

10" 

(5,50) 

(3.46 ±0.16) 

X 

10® 

(1.04 ±0.08) X 

10® 

(5.04 ± 

0.06) 

X 

10" 

(5.03 ± 

0.06) 

X 

10" 

(6,70) 

(1.01 ±0.30) 

X 

10" 

(2.21 ±0.03) X 

10" 

(5.12 ± 

0.15) 

X 

10" 

(5.05 ± 

0.22) 

X 

10" 

Ising model 

h = 3.01 













(2,20) 

(3.44 ± 1.85) 

X 

10^ 

5.97 ±0.09 


(1.29 ± 

0.13) 

X 

10" 

(1.29 ± 

0.14) 

X 

10" 

(3,30) 

(4.28 ±0.59) 

X 

10® 

(1.95 ±0.24) X 

10® 

(5.00 ± 

0.34) 

X 

10" 

(5.01 ± 

0.32) 

X 

10" 

(4,30) 

(1.28 ± 1.52) 

X 

10"® 

(5.81 ± 11.88) X 10"" 

(5.39 ± 

0.13) 

X 

10® 

(5.47 ± 

0.14) 

X 

10® 


Table I: Mean condition number of the norm matrix Nlr as well as of matrices L and R with their standard deviation. 
The numbers are for the case of positive approximant, hxing and without fixing the gauge, for the subtensor FFU. We show 
different bond dimensions [D, x) for iPEPS algorithm with time step S = 0.02. The mean is computed over ten time steps, 
for Heisenberg and Ising models. The transverse field for the Ising model is h = 3.01, close to criticality. 



Figure 10: (Color online) Local magnetization as a func¬ 
tion of the transverse held h in the quantum Ising model, with 
gauge (hlled circles) and without gauge (star symbols). The 
inset is a zoom around criticality. 


that, when the gauge-fixing is applied, both cost function 
and local fidelity tend to converge faster than the case 
without gauge fixing, see Fig. We also noted that 
speed-up in convergences becomes more significant as the 
bond dimension D of the iPEPS increases (not shown). 

We also compute the order parameter = (a^) for 
the quantum Ising model, and compare it for the cases 
with and without gauge fixing, see Fig. The results 
for both cases are quantitatively very similar deep in the 
gapped phases. However, when close to the quantum 
critical point, the results obtained with the iPEPS -I- 
gauge-fixing are better, and again we see that this ef¬ 
fect becomes more relevant for larger bond dimensions 
(e.g., we almost see no difference for D = 2, whereas 
for Z? = 3 we already see a clear difference). Be¬ 
sides, in Fig. we have plotted the correlation func¬ 
tion Szz{x) = at the critical point 

h = 3.044. We see that the simulation with gauge-fixing 
captures the correlation of the system better for the case 



X 


Figure 11: (Color online) Two-point correlator Szzix) of the 
quantum Ising model close to the critical point h = 3.044, 
computed with gauge (filled circles) and without gauge (star 
symbols). 


of bond dimensions {D,x) = (3,30). 

We have also applied the “improved” iPEPS algorithm 
to study the ground state of the Heisenberg model for dif¬ 
ferent bond dimensions {D,x)- For small bond dimen¬ 
sions D, we do not see much difference for the results 
with and without gauge fixing. But for large bond dimen¬ 
sions we observe that the ground state energy obtained 
using gauge fixing is better, see Fig. l^a). To quantify 
the overall error, we compare to the best result obtained 
from quantum Monte CarlcPH, eg = —0.669437(5). In our 
case, with gauge-fixing we obtain eg = —0.669309(2) and 
without gauge e„g = —0.669243(1), for {D,x) = (7,70). 
Besides, we also compute the staggered magnetization m 
as a function of the bond dimension, shown in Fig. lilb). 
Again, for large bond dimensions the calculations with 
the gauge-fixing become better than without the gauge. 
Our best values were obtained for {D,x) = (7,70), and 
are rUg = 0.33490 with gauge m„g = 0.33662 without 
gauge. This is to be compared with the Monte Carlo 
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D D 


Figure 12: (Color online) (a) Relative error of the energy per 
link e = (eo — e)/eo for the Heisenberg model, as a function of 
the iPEPS bond dimension D. Here eo = —0.669437(5) is the 
quantum Monte Carlo result.^ (b) Staggered magnetization 
m of the Heisenberg model as a function of the bond dimen¬ 
sion D, compared to the Monte Carlo result mo = 0.30703.® 
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Figure 13: (Color online) Total running time (in seconds) for 
five time steps in imaginary-time evolution with the Heisen¬ 
berg Hamiltonian, with the FU and FFU (-F gauge-Hxing) 
updates. The speed-up factor of the new update scheme with 
respect to the old one is shown in the inset. 


result. Too = 0.30703.1^ 

Finally, even if not mentioned explicitly at every calcu¬ 
lation, we remind that our results are obtained by using 
the FFU which is remarkably faster than the costly FU 
update, because the environment is not recomputed from 
scratch at every step. In Fig. [T^ we compare the actual 
running times for the FU and the FFU (-F gauge-fixing) 
schemes. The speed-up factor increases with increasing 
H, since more CTM steps are required to reach conver¬ 
gence of the environment tensors if the state is more en¬ 
tangled. In the present example we obtained a speed-up 
factor of up to ~ 30. This factor will be even larger in 


more strongly entangled systems (e.g. in fermionic sys¬ 
tems). In combination with the gauge-fixing the FFU 
approach is even better: we have seen that the overall 
improved algorithm is remarkably stable, as well as sub¬ 
stantially faster than the old version. 

One more comment is in order: one could naively ex¬ 
pect that after several updates with the FFU scheme the 
environment will have drifted away, so that it may need 
to be reinitialized by fully converging it. However, we 
find that the new update scheme is self correcting for the 
considered values of i5 and for the studied models (indeed, 
if 6 is small enough, the changes in the tensors are also 
expected to be small). One may expect, however, that 
for larger S the method is no longer self-correcting. This 
is a possibility that needs to be taken into account when 
implementing the algorithm in practice. But in any case, 
6 decreases throughout the evolution of the algorithm, 
so the self-correction happens naturally and in practice 
we never need to restart. For the models analyzed in 
our paper we never encounter such situation, but we can 
imagine that for more complex systems one may need 
to restart the environment from time to time in order 
to improve the results. This may be also an important 
difference between imaginary- and real-time evolutions. 
For imaginary evolutions, the algorithm is naturally self- 
correcting, since in the limit 5 —> 0 the actual time steps 
behave, in practice, as convergence steps for the environ¬ 
ment. However, this property may be lost for real-time 
evolutions if 6 is not small enough, so that several itera¬ 
tions may be needed at every step. But in any case, for 
real-time evolutions it is also important to recycle envi¬ 
ronment tensors at every iteration in order to save time, 
as well as to take care of a correct matching of all bond 
indices as done explicitly in the FFU. 


V. CONCLUSIONS 


In this paper we have presented how to improve the 
stability and efficiency of the iPEPS algorithm. We have 
discussed two improvements, namely the fast full update 
(FFU) and the gauge-fixing in the ALS sweep. In the 
FFU scheme the tensors in the effective environment are 
updated at every step, while keeping them compatible 
with the updated iPEPS tensors. This implies a speed¬ 
up by a large factor (up to ~ 30 in the present examples, 
and even larger in more strongly entangled systems) com¬ 
pared to the previous EU approach, where the environ¬ 
ment tensors have been recomputed from scratch at each 
time step. The gauge-fixing improves the conditioning in 
the ALS sweep at every update step, in the same way as 
was already shown for finite PEPS calculations ,1^ leading 
to a better stability and faster convergence. 

We have benchmarked the improved iPEPS algorithm 
with calculations for the quantum Ising and Heisenberg 
models on an infinite 2d square lattice, where we have 
seen that similar or slightly better accuracies can be ob¬ 
tained, substantially faster, and with more stable evo- 
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lutions, when compared to previous iPEPS calculations. 
This is particularly true for large bond dimensions and 
in highly-entangled systems, such as in the vicinity of 
quantum critical points. 

Technically, we have demonstrated the improved 
iPEPS algorithm for systems on an infinite square lattice 
and a 2-site unit cell, but the extension to other 2d lat¬ 
tices and bigger unit cells is straightforward. The method 
can also be extended to Hamiltonians with longer-range 
interactions. We expect that these improvements will be 
a significant step forward towards powerful tensor net¬ 
work calculations for challenging 2d systems. 
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